#!/usr/bin/env python
"""

    get_ncbi_tax_tree.py
    [--log_file PATH]
    [--verbose]

"""

################################################################################
#
#   get_ncbi_tax_tree
#
#
#   Copyright (c) 11/24/2009 Leo Goodstadt
#
#   Permission is hereby granted, free of charge, to any person obtaining a copy
#   of this software and associated documentation files (the "Software"), to deal
#   in the Software without restriction, including without limitation the rights
#   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#   copies of the Software, and to permit persons to whom the Software is
#   furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
#   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
#   THE SOFTWARE.
#################################################################################



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   options


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

import sys, os
import os.path


# add self to search path for testing
if __name__ == '__main__':
    exe_path = os.path.split(os.path.abspath(sys.argv[0]))[0]
    sys.path.append(os.path.abspath(os.path.join(exe_path,"..", "python_modules")))
    #sys.path.insert(0, "/net/cpp-group/Leo/inprogress/oss_projects/ruffus/installation/src")
    module_name = os.path.split(sys.argv[0])[1]
    module_name = os.path.splitext(module_name)[0];
else:
    module_name = __name__


if __name__ == '__main__':
    from optparse import OptionParser
    import StringIO

    parser = OptionParser(version="%prog 1.0", usage = "\n\n    %progs [options]")
    parser.add_option(  "-n", "--ncbi_tax_data_path", dest="ncbi_tax_data_path",
                        default = "/net/cpp-mirror/ncbi/taxonomy",
                        metavar="PATH",
                        type="string",
                        help="path to ncbi names.dmp and nodes.dmp files. Cached python data "
                             "will also be kept there.\n"
                             "names.dmp has fields of tax_id, name, "
                             "unique_name, name_class. Entries ignored unless name_class =="
                             "scientific name."
                             "nodes.dmp has tax_id, parent_tax_id as the first two fields.")
    parser.add_option(  "-s", "--species_file", dest="species_file",
                        default = "/net/cpp-mirror/ncbi/taxonomy/interesting.species",
                        metavar="FILE",
                        type="string",
                        help="name and path of file holding a list of species we are interested "
                            "in. Leave blank for the tree of life. ")

    #
    #   general options: verbosity / logging
    #
    parser.add_option("-v", "--verbose", dest = "verbose",
                      action="count", default=0,
                      help="Print more verbose messages for each additional verbose level.")
    parser.add_option("-L", "--log_file", dest="log_file",
                      metavar="FILE",
                      type="string",
                      help="Name and path of log file")
    parser.add_option("--skip_parameter_logging", dest="skip_parameter_logging",
                        action="store_true", default=False,
                        help="Do not print program parameters to log.")
    parser.add_option("-D", "--debug", dest="debug",
                        action="count", default=0,
                        help="Set default program parameters in debugging mode.")


    parser.set_defaults(log_file       = os.path.join(exe_path, "parse_ncbi.log"))


    # get help string
    f =StringIO.StringIO()
    parser.print_help(f)
    helpstr = f.getvalue()
    (options, remaining_args) = parser.parse_args()


    #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
    #                                             #
    #   Debug: Change these                       #
    #                                             #
    #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    if options.debug:
        options.verbose                 = 5
    #vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
    #                                             #
    #   Debug: Change these                       #
    #                                             #
    #^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

    #
    #   mandatory options
    #
    from options import check_mandatory_options
    mandatory_parameters = ["ncbi_tax_data_path"]
    check_mandatory_options (options, mandatory_parameters, helpstr)


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   imports


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#from json import dumps
from collections import defaultdict
import cPickle
import marshal
caching_method = marshal
#caching_method = cPickle



#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Functions


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Logger


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

if __name__ == '__main__':
    import logging
    from lg_program_logging import  setup_std_logging, MESSAGE
    from options import get_option_strings


    logger = logging.getLogger(module_name)
    setup_std_logging(logger, options.log_file, options.verbose)


    #
    #   log programme parameters
    #
    if not options.skip_parameter_logging:
        programme_name = os.path.split(sys.argv[0])[1]
        logger.info("%s %s" % (programme_name, " ".join(get_option_strings(parser, options))))

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Globals


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
tax_id_to_scientific_name       = dict()
scientific_name_to_tax_id       = dict()
tax_id_to_common_name           = dict()
scientific_name_to_common_name  = dict()
child_to_parent_tax_id          = dict()
parent_to_child_tax_id          = dict()
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Functions


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
import re
field_separator = re.compile("\t\|\t?");




#_________________________________________________________________________________________

#   _read_tax_names

#       get tax_id and taxonomy namess

#_________________________________________________________________________________________
def link_scientific_common_names (tax_id_to_common_name, tax_id_to_scientific_name):
    """
    Lookup table between common and scientific names
        use common name of parent if we don't have one
    """
    def link_names (scientific_name_to_common_name, tax_id1, tax_id2):
        scientific_name = tax_id_to_scientific_name[tax_id1]
        common_name     = tax_id_to_common_name[tax_id2].capitalize()
        scientific_name_to_common_name[scientific_name] = common_name
        ##DEBUGG
        #if tax_id == 9887 or tax_id == 9888:
        #    sys.stderr.write("Link parent names %s => %s\n" % (scientific_name, common_name))

    scientific_name_to_common_name = dict()
    for tax_id in tax_id_to_scientific_name:
        if tax_id in tax_id_to_common_name:
            link_names(scientific_name_to_common_name, tax_id, tax_id)
        #
        #   No common name, see if parent has one
        #
        elif tax_id in child_to_parent_tax_id:
            parent_tax_id = child_to_parent_tax_id[tax_id]
            if parent_tax_id in tax_id_to_common_name:
                link_names(scientific_name_to_common_name, tax_id, parent_tax_id)

    return scientific_name_to_common_name


def _read_tax_names (names_file, child_to_parent_tax_id):
    """
    Read Taxonomy names from ncbi names.dmp

    """



    #
    #
    if __name__ == '__main__':
        logger.info("Reading taxonomy names...")
    for line in open(names_file):
        fields = re.split(field_separator, line)
        if len(fields) < 4 or fields[0][0] == '#':
            continue
        tax_id, name_txt, unique_name, name_class  = fields[0:4]
        name_txt = name_txt.lower()
        tax_id = int(tax_id)
        if name_class == "scientific name":
            scientific_name = name_txt.capitalize()
            tax_id_to_scientific_name[tax_id] = scientific_name
            scientific_name_to_tax_id[scientific_name] = tax_id

            ##DEBUGG
            #if tax_id == 9887 or tax_id == 9888:
            #    sys.stderr.write("%s = %d\n" % (scientific_name, tax_id))

        # genbank common name takes precedence over "ordinary common names"
        elif name_class == "genbank common name":
            genbank_common_name = name_txt

            ##DEBUGG
            #if tax_id == 9887 or tax_id == 9888:
            #    if tax_id in tax_id_to_common_name:
            #        prev_name = tax_id_to_common_name[tax_id]
            #    else:
            #        prev_name = ""
            #    sys.stderr.write("genbank_common name for %d = %s => %s\n" % (tax_id, prev_name, common_name))

            #ignore if common name already superset of genbank common name
            if  (tax_id in tax_id_to_common_name and
                genbank_common_name in tax_id_to_common_name[tax_id]):

                ##DEBUGG
                #if tax_id == 9887 or tax_id == 9888:
                #    sys.stderr.write("Ignore because subset\n")

                continue

            tax_id_to_common_name[tax_id] = genbank_common_name

            ##DEBUGG
            #if tax_id == 9887 or tax_id == 9888:
            #    sys.stderr.write("genbank_common name for %d = %s\n" % (tax_id, tax_id_to_common_name[tax_id]))

        # only take first one: may be more than one
        # use longest superset of a name, even if that overwrites a genbank common name
        elif name_class == "common name":

            #DEBUGG
            #if tax_id == 9887 or tax_id == 9888:
            #    sys.stderr.write("common name for %d replaces? %d\n" % (tax_id, tax_id in tax_id_to_common_name))
            #
            if (not tax_id in tax_id_to_common_name or
                tax_id_to_common_name[tax_id] in name_txt) :

                ##DEBUGG
                #if tax_id == 9887 or tax_id == 9888:
                #    if tax_id in tax_id_to_common_name:
                #        sys.stderr.write("common name for %d %s => %s\n" % (tax_id, tax_id_to_common_name[tax_id], common_name))
                #    else:
                #        sys.stderr.write("common name for %d => %s\n" % (tax_id, common_name))

                tax_id_to_common_name[tax_id] = name_txt


    scientific_name_to_common_name = link_scientific_common_names (tax_id_to_common_name, tax_id_to_scientific_name)
    return (tax_id_to_scientific_name, tax_id_to_common_name, scientific_name_to_tax_id,
            scientific_name_to_common_name)





#_________________________________________________________________________________________

#   _read_tax_id_hierarchy

#       get tax_id and parent_tax_id

#_________________________________________________________________________________________
def _read_tax_id_hierarchy (nodes_file):
    """
    get tax_id and parent_tax_id

    """
    if __name__ == '__main__':
        logger.info("Reading node structure for taxonomy tree...")

    tmp_parent_to_child_tax_id = defaultdict(list)
    for line in open(nodes_file):
        fields = re.split(field_separator, line)
        if len(fields) < 3 or fields[0][0] == '#':
            continue
        if fields[0] == fields[1]:
            continue
        tmp_parent_to_child_tax_id[int(fields[1])].append(int(fields[0]))
        child_to_parent_tax_id[int(fields[0])] = int(fields[1])

    parent_to_child_tax_id.update(tmp_parent_to_child_tax_id)
    return parent_to_child_tax_id, child_to_parent_tax_id






#_________________________________________________________________________________________

#   _read_species_of_interest

#_________________________________________________________________________________________
def _read_species_of_interest (species_file):
    """
    Reading species whose hierarchy will be retained

    """
    logger.info("Reading species of interest...")

    species_of_interest = set()
    for line in open(species_file):
        line = line.rstrip()
        if not len(line) or line [0] == '#':
            continue
        species_of_interest.add(line)

    return species_of_interest







#_________________________________________________________________________________________

#   _restrict_data_to_hierarchy_of_interesting_species

#_________________________________________________________________________________________
def _restrict_data_to_hierarchy_of_interesting_species (   species_of_interest,
                                                            scientific_name_to_tax_id,
                                                            child_to_parent_tax_id,
                                                            parent_to_child_tax_id,
                                                            die_if_missing_species):
    """

    Mark all ancestors of interesting species for retention

    Mark all interesting species as terminal branches (deleting sub species if necessary)

    Will throw exception if any of the species are not in tax tree

    """
    tax_ids_in_hierarchy = set()

    for species in species_of_interest:
        if not species in scientific_name_to_tax_id:
            if die_if_missing_species:
                die_error("Missing species of interest [%s] in the NCBI taxonomy!" % species)
            else:
                continue

        tax_id = scientific_name_to_tax_id[species]
        # mark all ancestral tax_ids for retention
        retained_tax_id = tax_id
        tax_ids_in_hierarchy.add(retained_tax_id)
        while retained_tax_id in child_to_parent_tax_id:
            retained_tax_id = child_to_parent_tax_id[retained_tax_id]
            tax_ids_in_hierarchy.add(retained_tax_id)

        if tax_id in parent_to_child_tax_id:
            del parent_to_child_tax_id[tax_id]

    return tax_ids_in_hierarchy




#_________________________________________________________________________________________

#   read_ncbi_taxa_data

#_________________________________________________________________________________________
def _last_modified_date (filename):
    """
    retrieve last modified date of file
    """
    stats = os.stat(filename)
    return stats[8]

def read_ncbi_taxa_data(force_reread = False):
    """
    Read ncbi taxa either from previously cached or by parsing files

        sets variables for the whole module.
        Just call this function again after data files changed

    """

    global tax_id_to_scientific_name
    global tax_id_to_common_name
    global scientific_name_to_tax_id
    global parent_to_child_tax_id
    global child_to_parent_tax_id
    global scientific_name_to_common_name



    CACHED_DATA = os.path.join(NCBI_TAXA_DIR, "cached_ncbi_taxa")

    #
    #   cached data exists and is younger than the names or nodes files
    #
    if  (not force_reread and
         os.path.exists(CACHED_DATA) and
         _last_modified_date(CACHED_DATA) > _last_modified_date(NAMES_FILE_NAME) and
         _last_modified_date(CACHED_DATA) > _last_modified_date(NODES_FILE_NAME)):

        if __name__ == '__main__':
            logger.info("Get cached ncbi taxa data from \n[%s]" % CACHED_DATA)
        cache_file = open(CACHED_DATA, 'rb')
        parent_to_child_tax_id          = caching_method.load(cache_file)
        child_to_parent_tax_id          = caching_method.load(cache_file)

        tax_id_to_scientific_name       = caching_method.load(cache_file)
        scientific_name_to_tax_id       = dict((v,k) for (k, v) in tax_id_to_scientific_name.iteritems())
        #scientific_name_to_tax_id      = caching_method.load(cache_file)
        tax_id_to_common_name           = caching_method.load(cache_file)

        scientific_name_to_common_name = link_scientific_common_names (tax_id_to_common_name, tax_id_to_scientific_name)
        #scientific_name_to_common_name  = caching_method.load(cache_file)
        cache_file.close()

    #
    #   save to cache
    #
    else:
        # parse files and read data
        (parent_to_child_tax_id,
         child_to_parent_tax_id)            = _read_tax_id_hierarchy(NODES_FILE_NAME)
        (tax_id_to_scientific_name, scientific_name_to_tax_id,
         scientific_name_to_common_name,
         tax_id_to_common_name)             = _read_tax_names (NAMES_FILE_NAME,
                                                                child_to_parent_tax_id)

        #   cache
        if __name__ == '__main__':
            logger.info("Cache ncbi taxa data to \n[%s]" % CACHED_DATA)
        cache_file = open(CACHED_DATA, 'wb')
        caching_method.dump(parent_to_child_tax_id,         cache_file, 1)
        caching_method.dump(child_to_parent_tax_id,         cache_file, 1)
        caching_method.dump(tax_id_to_scientific_name,      cache_file, 1)
        caching_method.dump(tax_id_to_common_name,          cache_file, 1)
        cache_file.close()

#_________________________________________________________________________________________

#   get_descendants

#_________________________________________________________________________________________
def get_descendants(ancestral_tax_id, tax_ids_of_interest=None):
    """
    Read ncbi taxa either from previously cached or by parsing files

        sets variables for the whole module.
        Just call this function again after data files changed

    """
    if tax_ids_of_interest == None:
        tax_ids_of_interest = tax_id_to_scientific_name.keys()

    descendants = set([ancestral_tax_id])
    non_descendants = set()

    for tax_id in tax_ids_of_interest:


        if tax_id == ancestral_tax_id:
            descendants.add(tax_id)
            continue

        curr_tax_id = tax_id
        tax_ids_in_play = [curr_tax_id]

        while curr_tax_id in child_to_parent_tax_id :
            if curr_tax_id in descendants:
                descendants.update(tax_ids_in_play)
                tax_ids_in_play = []
                break
            elif curr_tax_id in non_descendants:
                break
            # keep going
            else:
                curr_tax_id = child_to_parent_tax_id[curr_tax_id]
                tax_ids_in_play.append(curr_tax_id)

        # got to the top of tree and not found ancestral so must be non-descendants
        non_descendants.update(tax_ids_in_play)

    # return set which is common to both
    return descendants & set(tax_ids_of_interest)


def get_placental_species(species_of_interest=None):
    """
    Read ncbi taxa either from previously cached or by parsing files

        sets variables for the whole module.
        Just call this function again after data files changed

    """
    ancestral_placental_tax_id = scientific_name_to_tax_id["Eutheria"]

    #
    #   convert species of interest to tax_ids
    #
    if len(species_of_interest):
        tax_ids_of_interest = [scientific_name_to_tax_id[n] for n in species_of_interest]
    else:
        tax_ids_of_interest = set(tax_id_to_scientific_name.keys())
    all_placental_tax_ids = get_descendants(ancestral_placental_tax_id,
                                            tax_ids_of_interest)
    return sorted(tax_id_to_scientific_name[t] for t in all_placental_tax_ids)

def get_taxa_of_interest(taxa_of_interest, species_of_interest=None):
    """                                                                                                                                                                                                          
    Read ncbi taxa either from previously cached or by parsing files                                                                                                                                             
                                                                                                                                                                                                                 
        sets variables for the whole module.                                                                                                                                                                     
        Just call this function again after data files changed                                                                                                                                                   
                                                                                                                                                                                                                 
    """
    ancestral_tax_id = scientific_name_to_tax_id[taxa_of_interest]

    #                                                                                                                                                                                                            
    #   convert species of interest to tax_ids                                                                                                                                                                   
    #                                                             
    if len(species_of_interest):
        tax_ids_of_interest = [scientific_name_to_tax_id[n] for n in species_of_interest]
    else:
        tax_ids_of_interest = set(tax_id_to_scientific_name.keys())

    all_interesting_tax_ids = get_descendants(ancestral_tax_id,
                                            tax_ids_of_interest)
    return sorted(tax_id_to_scientific_name[t] for t in all_interesting_tax_ids)



from Tree import Tree
from Bio.Nexus import Nodes

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
#
#   Functions for adding branches and leaves to Nexus trees
#
def _add_branch (tree, parent_id):
    nd=tree.dataclass()
    sn=Nodes.Node(nd)
    tree.add(sn,parent_id)
    return sn.id

def _add_leaf (tree, leaf_name, parent_id):
    nd=tree.dataclass()
    nd.taxon=leaf_name
    leaf=Nodes.Node(nd)
    tree.add(leaf,parent_id)


#_________________________________________________________________________________________

#   _create_tree_by_recursion

#_________________________________________________________________________________________
def _create_tree_by_recursion (tree, parent_to_child_tax_id,
                        tax_id_to_scientific_name, tax_ids_in_hierarchy, parent, current_tax_id):

    children_tax_ids = parent_to_child_tax_id[current_tax_id]
    for child_tax_id in tax_ids_in_hierarchy.intersection(children_tax_ids):
        if not child_tax_id in parent_to_child_tax_id:
            if not tax_id_to_scientific_name[child_tax_id]:
                raise "missing %d" % child_tax_id
            _add_leaf(tree, tax_id_to_scientific_name[child_tax_id], parent)
        else:
            grandchild_tax_ids = tax_ids_in_hierarchy.intersection(parent_to_child_tax_id[child_tax_id])
            if len(grandchild_tax_ids) > 1:
                new_branch = _add_branch(tree, parent)
                _create_tree_by_recursion(tree, parent_to_child_tax_id, tax_id_to_scientific_name,
                                tax_ids_in_hierarchy, new_branch, child_tax_id)
            # skip child if grandchild has a single parent
            else:
                _create_tree_by_recursion(tree, parent_to_child_tax_id, tax_id_to_scientific_name,
                                tax_ids_in_hierarchy, parent, child_tax_id)

#
#
#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


if __name__ == '__main__':
    NCBI_TAXA_DIR = options.ncbi_tax_data_path
else:
    NCBI_TAXA_DIR = "/net/cpp-mirror/ncbi/taxonomy"


NODES_FILE_NAME = os.path.join(NCBI_TAXA_DIR, "nodes.dmp")
NAMES_FILE_NAME = os.path.join(NCBI_TAXA_DIR, "names.dmp")

def get_ncbi_tree (species_of_interest, exception_if_species_missing_from_ncbi = False):
    """
    This is the main callable
    Get taxonomy tree from ncbi either by parsing or by retrieving from cache
    Restrict the tree to the species of interest
    """

    if len(species_of_interest):
        tax_ids_in_hierarchy = _restrict_data_to_hierarchy_of_interesting_species (set(species_of_interest),
                                                                               scientific_name_to_tax_id,
                                                                               child_to_parent_tax_id,
                                                                               parent_to_child_tax_id,
                                                                               exception_if_species_missing_from_ncbi)
    else:
        tax_ids_in_hierarchy = set(tax_id_to_scientific_name.keys())

    species_tree = Tree()
    sys.setrecursionlimit(100)
    _create_tree_by_recursion (species_tree, parent_to_child_tax_id, tax_id_to_scientific_name,
                                tax_ids_in_hierarchy, species_tree.root, 1)
    return species_tree




if __name__ != '__main__':
    read_ncbi_taxa_data()
else:
    if not options.debug:
        read_ncbi_taxa_data()




#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Main logic


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
if __name__ == '__main__':
    if options.debug:
        import unittest
        class Test_get_ncbi_tax_tree(unittest.TestCase):

            #       self.assertEqual(self.seq, range(10))
            #       self.assert_(element in self.seq)
            #       self.assertRaises(ValueError, random.sample, self.seq, 20)



            def test_function(self):
                """
                    test
                """
                read_ncbi_taxa_data()
                #print "\n"
                #print "-"*80
                #for i, (c, s) in enumerate(common_name_to_scientific_name.iteritems()):
                #    print "%40s\t%40s" % (c,s)
                #    if i> 20:
                #        break
                #print "-"*40
                #for i, (c, s) in enumerate(scientific_name_to_common_name.iteritems()):
                #    print "%40s\t%40s" % (c,s)
                #    if i> 20:
                #        break
                #print "-"*40
                #print "\n"
                species_list = ["Microcebus murinus"         ,
                                "Callithrix jacchus"         ,
                                "Sus scrofa"                 ,
                                "Gorilla gorilla gorilla"    ,
                                "Macaca mulatta"             ,
                                "Ateles geoffroyi"           ,
                                "Homo sapiens"               ,
                                "Pan troglodytes"            ,
                                "Mus musculus"               ,
                                "Mustela putorius furo"      ,
                                "Cavia porcellus"            ,
                                "Cavia porcellus"            ,
                                "Taeniopygia guttata"        ,
                                "Gallus gallus"              ,
                                "Monodelphis domestica"      ,
                                "Ornithorhynchus anatinus"   ,
                                "Tachyglossus aculeatus"     ]

                species_tree = get_ncbi_tree (species_list, True)
                import TreeTools
                tree_str = TreeTools.Tree2Newick(species_tree, False, False);
                self.assertEqual(tree_str,
                                    "(((Gallus gallus,Taeniopygia guttata),((((Sus scrofa,Mustela putorius furo),((Mus musculus,Cavia porcellus),(((Macaca mulatta,(Gorilla gorilla gorilla,Pan troglodytes,Homo sapiens)),(Callithrix jacchus,Ateles geoffroyi)),Microcebus murinus))),Monodelphis domestica),(Ornithorhynchus anatinus,Tachyglossus aculeatus))));")

            def test_get_placental(self):
                species_list = ["Microcebus murinus"         ,
                                "Callithrix jacchus"         ,
                                "Sus scrofa"                 ,
                                "Gorilla gorilla gorilla"    ,
                                "Macaca mulatta"             ,
                                "Ateles geoffroyi"           ,
                                "Homo sapiens"               ,
                                "Pan troglodytes"            ,
                                "Mus musculus"               ,
                                "Mustela putorius furo"      ,
                                "Cavia porcellus"            ,
                                "Cavia porcellus"            ,
                                "Taeniopygia guttata"        ,
                                "Gallus gallus"              ,
                                "Monodelphis domestica"      ,
                                "Ornithorhynchus anatinus"   ,
                                "Tachyglossus aculeatus"     ]
                placental_species_list = get_placental_species(species_list)
                self.assertEqual(placental_species_list,
                                    ['Ateles geoffroyi',
                                    'Callithrix jacchus',
                                    'Cavia porcellus',
                                    'Gorilla gorilla gorilla',
                                    'Homo sapiens',
                                    'Macaca mulatta',
                                    'Microcebus murinus',
                                    'Mus musculus',
                                    'Mustela putorius furo',
                                    'Pan troglodytes',
                                    'Sus scrofa'])


        #
        #   call unit test without parameters
        #

        if sys.argv.count("--debug"):
            sys.argv.remove("--debug")
        unittest.main()


    else:
        if options.species_file:
            logger.info("Restricting data to hierarchy of interesting species...")
            species_of_interest = _read_species_of_interest (options.species_file)
            tax_ids_in_hierarchy = _restrict_data_to_hierarchy_of_interesting_species (species_of_interest,
                                                                                   scientific_name_to_tax_id,
                                                                                   child_to_parent_tax_id,
                                                                                   parent_to_child_tax_id,
                                                                                   True)
        else:
            tax_ids_in_hierarchy = set(tax_id_to_scientific_name.keys())
        logger.info("Making tree...")
        species_tree = Tree()
        sys.setrecursionlimit(100)
        _create_tree_by_recursion (species_tree, parent_to_child_tax_id, tax_id_to_scientific_name,
                                    tax_ids_in_hierarchy, species_tree.root, 1)
        import TreeTools
        print TreeTools.Tree2Newick(species_tree, False, False);





